# set working directory
setwd("~/GitHub/CorticalFoldingHumans_deaging_harmonization")
set.seed(1)
# load packages ----
library(tidyverse)
library(readxl)
library(stringr)
library(lubridate)
library(kableExtra)
library(ggpubr)
library(broom)
library(MASS)
library(irr)
library(effects)
library(lme4)
library(lsmeans)
# load datasets ----
# removing OASIS as it does not have lobes information
data <- read_csv("data.csv") %>%
filter(Sample != "OASIS")
# estimate localGI, k, K, S and I
# correct for gaussian curvature for lobes estimates
# estimate localGI, k, K, S and I for lobes
data <- data %>%
mutate(
logAvgThickness = log10(AvgThickness),
logTotalArea = log10(TotalArea),
logExposedArea = log10(ExposedArea),
localGI = TotalArea / ExposedArea,
k = sqrt(AvgThickness) * TotalArea / (ExposedArea ^ 1.25),
K = 1 / 4 * log10(AvgThickness^ 2) + log10(TotalArea) - 5 / 4 * log10(ExposedArea),
S = 3 / 2 * log10(TotalArea) + 3 / 4 * log10(ExposedArea) - 9 / 4 * log10(AvgThickness^2),
I = log10(TotalArea) + log10(ExposedArea) + log10(AvgThickness^ 2),
Knorm = K / sqrt(1 + (1 / 4) ^ 2 + (5 / 4) ^ 2),
Snorm = S / sqrt((3 / 2) ^ 2 + (3 / 4) ^ 2 + (9 / 4) ^ 2),
Inorm = I / sqrt(1 ^ 2 + 1 ^ 2 + 1 ^ 2),
c = as.double(ifelse(ROI == "hemisphere", NA, 4 * pi / GaussianCurvature)),
TotalArea_corrected = ifelse(ROI == "hemisphere", TotalArea, TotalArea * c),
ExposedArea_corrected = ifelse(ROI == "hemisphere", ExposedArea, ExposedArea * c),
logTotalArea_corrected = log10(TotalArea_corrected),
logExposedArea_corrected = log10(ExposedArea_corrected),
localGI_corrected = TotalArea_corrected / ExposedArea_corrected,
k_corrected = sqrt(AvgThickness) * TotalArea_corrected / (ExposedArea_corrected ^ 1.25),
K_corrected = log10(TotalArea_corrected) - 5 / 4 * log10(ExposedArea_corrected) + 1 / 4 * log10(AvgThickness^ 2),
I_corrected = log10(TotalArea_corrected) + log10(ExposedArea_corrected) + log10(AvgThickness^ 2),
S_corrected = 3 / 2 * log10(TotalArea_corrected) + 3 / 4 * log10(ExposedArea_corrected) - 9 / 4 * log10(AvgThickness^ 2)
)
# Define variables as factor ----
data$ROI <- as.factor(data$ROI)
data$Diagnostic <- as.factor(data$Diagnostic)
data$Sample <- as.factor(data$Sample)
data$Gender <- as.factor(data$Gender)
# select CTL subjects ----
data <- data %>% filter(Diagnostic == "CTL")
# target age for deaging ----
Age.cor = 0
| Sample | N | age | age_range | AvgT | AT | AE |
|---|---|---|---|---|---|---|
| ADNI | 868 | 75 ± 6.5 | 56 ; 96 | 2.4 ± 0.11 | 97000 ± 8400 | 39000 ± 2700 |
| AHEAD | 100 | 42 ± 19 | 24 ; 76 | 2.6 ± 0.15 | 110000 ± 10000 | 39000 ± 2700 |
| AOMICPIOP1 | 208 | 22 ± 1.8 | 18 ; 26 | 2.6 ± 0.084 | 110000 ± 10000 | 40000 ± 2700 |
| AOMICPIOP2 | 224 | 22 ± 1.8 | 18 ; 26 | 2.6 ± 0.089 | 110000 ± 9100 | 39000 ± 2700 |
| HCPr900 | 881 | 29 ± 3.6 | 24 ; 37 | 2.7 ± 0.09 | 110000 ± 11000 | 40000 ± 3000 |
| IXI | 563 | 49 ± 16 | 20 ; 86 | 2.5 ± 0.14 | 97000 ± 10000 | 39000 ± 3100 |
| NKI | 168 | 34 ± 19 | 4 ; 85 | 2.5 ± 0.14 | 110000 ± 11000 | 40000 ± 2800 |
| Sample | N | age | age_range | K | S | I |
|---|---|---|---|---|---|---|
| ADNI | 868 | 75 ± 6.5 | 56 ; 96 | -0.56 ± 0.017 | 9.2 ± 0.13 | 10 ± 0.074 |
| AHEAD | 100 | 42 ± 19 | 24 ; 76 | -0.5 ± 0.023 | 9.2 ± 0.13 | 10 ± 0.094 |
| AOMICPIOP1 | 208 | 22 ± 1.8 | 18 ; 26 | -0.49 ± 0.014 | 9.1 ± 0.094 | 10 ± 0.078 |
| AOMICPIOP2 | 224 | 22 ± 1.8 | 18 ; 26 | -0.51 ± 0.013 | 9.1 ± 0.093 | 10 ± 0.078 |
| HCPr900 | 881 | 29 ± 3.6 | 24 ; 37 | -0.5 ± 0.013 | 9.1 ± 0.11 | 10 ± 0.082 |
| IXI | 563 | 49 ± 16 | 20 ; 86 | -0.56 ± 0.022 | 9.2 ± 0.14 | 10 ± 0.1 |
| NKI | 168 | 34 ± 19 | 4 ; 85 | -0.52 ± 0.027 | 9.2 ± 0.12 | 10 ± 0.094 |
Visual verification of morphological variables behaviour with aging. Considering not all brain regions ages in the same rate, we split it in the desired ROI.
ap1_b <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = TotalArea)) +
geom_point(aes(alpha =0.4, color = Sample, fill = Sample ),
size = 4,
pch = 21,
colour = 'white') +
geom_smooth(method = "lm", color = "black") +
stat_cor() +
guides(alpha= "none") +
theme_pubr()+
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25)))
ap1_b
ap1_c <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = ExposedArea)) +
geom_point(aes(alpha =0.4, color = Sample, fill = Sample ),
size = 4,
pch = 21,
colour = 'white') +
geom_smooth(method = "lm", color = "black") +
stat_cor() +
guides(alpha= "none") +
theme_pubr()+
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25)))
ap1_c
ap1_e <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = S)) +
geom_point(aes(alpha =0.4, color = Sample, fill = Sample ),
size = 4,
pch = 21,
colour = 'white') +
geom_smooth(method = "lm", color = "black") +
stat_cor() +
guides(alpha= "none") +
theme_pubr()+
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25)))
ap1_e
ap1_f <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = I)) +
geom_point(aes(alpha =0.4, color = Sample, fill = Sample),
size = 4,
pch = 21,
colour = 'white') +
geom_smooth(method = "lm", color = "black") +
stat_cor() +
guides(alpha= "none") +
theme_pubr()+
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25)))
ap1_f
# Generate and save figures
ap1 <-
ggarrange(ap1_a, ap1_b, ap1_c, ap1_d, ap1_e, ap1_f, labels = c("A","B", "C", "D", "E", "F"), nrow = 3, ncol = 2, font.label = list(size = 11), common.legend = TRUE, legend = "top")
ggsave("ap1.pdf", ap1, dpi=1200, width = 17.8, height = 17.8, units = "cm", device = "pdf")
ggplot(filter(data, ROI != "hemisphere"), aes(x = Age, y = TotalArea)) +
geom_point(aes(alpha =0.4, color = Sample, fill = Sample),
size = 4,
pch = 21,
colour = 'white') +
geom_smooth(method = "lm", color = "black") +
stat_cor() +
facet_grid(ROI ~ .) +
guides(alpha= "none") +
theme_pubr()+
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25)))
ggplot(filter(data, ROI != "hemisphere"), aes(x = Age, y = ExposedArea)) +
geom_point(aes(alpha =0.4, color = Sample, fill = Sample),
size = 4,
pch = 21,
colour = 'white') +
geom_smooth(method = "lm", color = "black") +
stat_cor() +
facet_grid(ROI ~ .) +
guides(alpha= "none") +
theme_pubr()+
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25)))
ggplot(filter(data, ROI != "hemisphere"), aes(x = Age, y = S)) +
geom_point(aes(alpha =0.4, color = Sample, fill = Sample),
size = 4,
pch = 21,
colour = 'white') +
geom_smooth(method = "lm", color = "black") +
stat_cor() +
facet_grid(ROI ~ .) +
guides(alpha= "none") +
theme_pubr() +
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25)))
ggplot(filter(data, ROI != "hemisphere"), aes(x = Age, y = I)) +
geom_point(aes(alpha =0.4, color = Sample, fill = Sample),
size = 4,
pch = 21,
colour = 'white') +
geom_smooth(method = "lm", color = "black") +
stat_cor() +
facet_grid(ROI ~ .) +
guides(alpha= "none") +
theme_pubr() +
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25)))
We remove the age effect at the base measurements, Total Area, Exposed Area and Cortical Thickness.
For the deaging process, we use Robust Linear Regressions, as it give weights to each data point and reduce the effects of outliers.
The deaging process NEED to be based on CONTROL values. Age also affects brain structure, as seens in the figures above. If the deaging process was based on all subjects, we could remove the disease effect across aging. Also, the deaging process is unique for each ROI. Therefore, we regress the variable for each ROI in an independent process here.
Further, we also need to separate for each Sample. All humans tends to have the same behaviour, but heterogeneous acquisition and processing could imply in a Type II Error addition in the results. This will be handle by the harmonization process.
## AvgThickness ---
# select only CTL subjects with logAvgThickness values.
# divide by each ROI and sample
# perform the robust linear regression and visualize all results as a data frame
decay_AvgThickness <-
filter(data, Diagnostic == "CTL", !is.na(AvgThickness), !is.nan(AvgThickness), !is.infinite(AvgThickness)) %>%
droplevels() %>%
group_by(ROI, Sample) %>%
do(fit_decay_AvgThickness = tidy(rlm(AvgThickness ~ Age, data = .), conf.int=TRUE)) %>% unnest(cols = c(fit_decay_AvgThickness))
# display the results for example
decay_AvgThickness %>%
filter(ROI == "hemisphere") %>%
kable() %>%
kable_styling()
| ROI | Sample | term | estimate | std.error | statistic | conf.low | conf.high |
|---|---|---|---|---|---|---|---|
| hemisphere | ADNI | (Intercept) | 2.8705639 | 0.0297791 | 96.395401 | 2.8121980 | 2.9289298 |
| hemisphere | ADNI | Age | -0.0068515 | 0.0003969 | -17.261100 | -0.0076294 | -0.0060735 |
| hemisphere | AHEAD | (Intercept) | 2.7944993 | 0.0193034 | 144.767398 | 2.7566654 | 2.8323332 |
| hemisphere | AHEAD | Age | -0.0057785 | 0.0004165 | -13.874348 | -0.0065948 | -0.0049622 |
| hemisphere | AOMICPIOP1 | (Intercept) | 2.7919001 | 0.0533391 | 52.342438 | 2.6873573 | 2.8964429 |
| hemisphere | AOMICPIOP1 | Age | -0.0076430 | 0.0023971 | -3.188472 | -0.0123412 | -0.0029448 |
| hemisphere | AOMICPIOP2 | (Intercept) | 2.9997519 | 0.0518957 | 57.803467 | 2.8980382 | 3.1014656 |
| hemisphere | AOMICPIOP2 | Age | -0.0177330 | 0.0023559 | -7.527025 | -0.0223506 | -0.0131155 |
| hemisphere | HCPr900 | (Intercept) | 2.7921977 | 0.0166452 | 167.747567 | 2.7595736 | 2.8248218 |
| hemisphere | HCPr900 | Age | -0.0035976 | 0.0005716 | -6.294081 | -0.0047179 | -0.0024773 |
| hemisphere | IXI | (Intercept) | 2.7075605 | 0.0095766 | 282.725861 | 2.6887907 | 2.7263303 |
| hemisphere | IXI | Age | -0.0048540 | 0.0001865 | -26.032197 | -0.0052194 | -0.0044885 |
| hemisphere | NKI | (Intercept) | 2.6985922 | 0.0122927 | 219.528357 | 2.6744990 | 2.7226854 |
| hemisphere | NKI | Age | -0.0045080 | 0.0003168 | -14.230010 | -0.0051289 | -0.0038871 |
For removing the age effect, we are insterested in the slope. The slope represents the, supposed for this task, constant rate of decreasing of each variable. Future works can improve this process by describing this behaviour with higher levels functions. We will then create a new variable, called c, as the correction constant. To describe fully, we also include the standard error of each regression slope.
decay_AvgThickness <- filter(decay_AvgThickness,term == "Age") %>%
mutate(c_AvgThickness = estimate,
std_error_c_AvgThickness = std.error) %>%
dplyr::select(c(ROI, Sample, c_AvgThickness, std_error_c_AvgThickness))
This process will be repeated for the Total Area and Exposed Area.
## TotalArea ---
decay_TotalArea <- filter(data, Diagnostic == "CTL", !is.na(TotalArea), !is.nan(TotalArea), !is.infinite(TotalArea)) %>%
droplevels() %>%
group_by(ROI, Sample) %>%
do(fit_decay_TotalArea = tidy(rlm(TotalArea ~ Age, data = .),conf.int=TRUE)) %>%
unnest(cols = c(fit_decay_TotalArea))
decay_TotalArea <- filter(decay_TotalArea, term == "Age") %>%
mutate(c_TotalArea = estimate,
std_error_c_TotalArea = std.error) %>%
dplyr::select(c(ROI, Sample, c_TotalArea, std_error_c_TotalArea))
## ExposedArea ---
decay_ExposedArea <- filter(data, Diagnostic == "CTL", !is.na(ExposedArea), !is.nan(ExposedArea), !is.infinite(ExposedArea)) %>%
droplevels() %>%
group_by(ROI, Sample) %>%
do(fit_decay_ExposedArea = tidy(rlm(ExposedArea ~ Age, data = .), conf.int = TRUE)) %>%
unnest(cols = c(fit_decay_ExposedArea))
decay_ExposedArea <- filter(decay_ExposedArea, term == "Age") %>%
mutate(c_ExposedArea = estimate,
std_error_c_ExposedArea = std.error) %>%
dplyr::select(c(ROI, Sample, c_ExposedArea, std_error_c_ExposedArea))
With the correction constants, we will correct every variable and then estimate K, S and I and its normalized values. It is possible to specify any Age reference for the deaging.
Considering the linear function \(y = ax + b\), we will estimate the new y trying to have an slope closer to zero (as the original slope will be subtract by its own approximation extracted from the robust linear regression). We have them \(y\_{new} = y - a(x - x\_{new})\), in our terms, \(y\_{new} = y - c(x - Age.cor)\).
data <- full_join(data, decay_AvgThickness) %>%
full_join(decay_TotalArea) %>%
full_join(decay_ExposedArea) %>%
mutate(
AvgThickness_age_decay = AvgThickness - c_AvgThickness * (Age - Age.cor),
TotalArea_age_decay = TotalArea - c_TotalArea * (Age - Age.cor),
ExposedArea_age_decay = ExposedArea - c_ExposedArea * (Age - Age.cor),
K_age_decay = log10(TotalArea_age_decay) + 1/4*log10(AvgThickness_age_decay^2) - 5/4*log10(ExposedArea_age_decay),
I_age_decay = log10(TotalArea_age_decay) + log10(ExposedArea_age_decay) + log10(AvgThickness_age_decay^2),
S_age_decay = 3/2*log10(TotalArea_age_decay) + 3/4*log10(ExposedArea_age_decay) - 9/4*log10(AvgThickness_age_decay^2),
Knorm_age_decay = K_age_decay/sqrt(1 + (1/4)^2 + (5/4)^2),
Snorm_age_decay = S_age_decay/sqrt((3/2)^2 + (3/4)^2 + (9/4)^2),
Inorm_age_decay = I_age_decay/sqrt(1^2+1^2+1^2))
To compare the values, we will have the raw measurements in light gray, and the new values in dark blue
Let’s change to 25 years old. At this age we expect the brain already finished development and did not started aging.
Age.cor = 25 #
data <- data %>%
mutate(
AvgThickness_age_decay = AvgThickness - c_AvgThickness * (Age - Age.cor),
TotalArea_age_decay = TotalArea - c_TotalArea * (Age - Age.cor),
ExposedArea_age_decay = ExposedArea - c_ExposedArea * (Age - Age.cor),
K_age_decay = log10(TotalArea_age_decay) + 1/4*log10(AvgThickness_age_decay^2) - 5/4*log10(ExposedArea_age_decay),
I_age_decay = log10(TotalArea_age_decay) + log10(ExposedArea_age_decay) + log10(AvgThickness_age_decay^2),
S_age_decay = 3/2*log10(TotalArea_age_decay) + 3/4*log10(ExposedArea_age_decay) - 9/4*log10(AvgThickness_age_decay^2),
Knorm_age_decay = K_age_decay/sqrt(1 + (1/4)^2 + (-5/4)^2),
Snorm_age_decay = S_age_decay/sqrt((3/2)^2 + (3/4)^2 + (9/4)^2),
Inorm_age_decay = I_age_decay/sqrt(1^2+1^2+1^2))
How changing the correction Age affects of Cortical Thickness and K deaging:
The harmonization procedure is useful to combine multiple datasets to improve the statistical power, expand the findings to other Samples and studies and validating the findings of a local Sample to global results. The heterogeneous methodology limits these comparisons, so we will handle the differences to work the data.
The harmonization procedure follow the same premises as the deaging, but here we assume all samples follow the same trend. We will them remove the Type II Error by subtracting for each datapoint the residue of its sample.
If we extrapolate the regression lines we have:
We can see how the intercepts at 0 years old are different for the three Samples. With the harmonization procedure, our goal is to have the same intercept for all Samples, because we are assuming this difference is due to the different acquisition and processing materials and methods.
To do this we need to consider the Sample as a random effect in a Linear Mixed Model. The residue will be then the difference between each Sample intercept and a general intercept for all datapoints (black dashed line).
Again, our reference are the Control subjects. To keep it simple, we will harmonize the samples only for the hemisphere in this example. The code including the lobes will be available right after.
data_rate <-
filter(
data,
ROI == "hemisphere",
Diagnostic == "CTL"
)
data_rate$Sample <- as.factor(data_rate$Sample)
As in the deaging, we harmonize the raw measurements. We will use the Cortical Thickness as example again.
## AvgThickness \~ Age ----
m.1 <- lme4::lmer(AvgThickness ~ Age + (1|Sample), data = data_rate) # run the linear mixed model
as_tibble(ranef(m.1)) # random effects in tibble format
| grpvar | term | grp | condval | condsd |
|---|---|---|---|---|
| Sample | (Intercept) | ADNI | 0.0008163 | 0.0024061 |
| Sample | (Intercept) | AHEAD | 0.0210192 | 0.0070004 |
| Sample | (Intercept) | AOMICPIOP1 | -0.0149690 | 0.0048896 |
| Sample | (Intercept) | AOMICPIOP2 | -0.0289463 | 0.0047141 |
| Sample | (Intercept) | HCPr900 | 0.0833446 | 0.0023884 |
| Sample | (Intercept) | IXI | -0.0295801 | 0.0029849 |
| Sample | (Intercept) | NKI | -0.0316847 | 0.0054318 |
ap5_a <- ggplot(as_tibble(ranef(m.1)), aes(x = grp, y = condval)) +
geom_pointrange(aes(ymin = condval - condsd, ymax = condval + condsd)) +
geom_hline(yintercept = 0, linetype = "dashed") +
theme_pubr() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(x = "Sample", y = "Residual (standard deviation)")
ap5_a
ggsave("ap5_a.pdf", ap5_a, dpi=1200, width = 8.7, height = 8.7, units = "cm", device = "pdf")
re <- as_tibble(ranef(m.1)) %>% # extract the random effects
mutate(T_shift = condval, Sample = grp) %>% # the condval is the intercept difference for each sample. it is also a constant. here we will refer to it as shift
dplyr::select(-c(condval, grpvar, term, condsd, grp))
data <- full_join(data, re) %>% mutate(AvgThickness_shiftc = AvgThickness - T_shift) # includign the constant in the dataset and estimanting the harmonizaed values
Condval (in the figure called “Residual”) equal to zero means the sample intercept is the same as the general regressing line.
## TotalArea \~ Age ----
m.1 <- lme4::lmer(TotalArea ~ Age + (1|Sample), data = data_rate)
re <- as_tibble(ranef(m.1)) %>%
filter(grpvar == "Sample") %>%
mutate(AT_shift = condval, Sample = grp) %>%
dplyr::select(-c(condval, grpvar, term, condsd, grp))
data <- full_join(data, re) %>% mutate(TotalArea_shiftc = TotalArea - AT_shift)
## ExposedArea \~ Age ----
m.1 <- lme4::lmer(ExposedArea ~ Age + (1|Sample), data = data_rate)
re <- as_tibble(ranef(m.1)) %>%
filter(grpvar == "Sample") %>%
mutate(AE_shift = condval, Sample = grp) %>%
dplyr::select(-c(condval, grpvar, term, condsd, grp))
data <- full_join(data, re) %>% mutate(ExposedArea_shiftc = ExposedArea - AE_shift)
Re-estimating values:
data <- data %>%
mutate(
AvgThickness_shiftc = AvgThickness - c_AvgThickness * (Age - Age.cor),
TotalArea_shiftc = TotalArea - c_TotalArea * (Age - Age.cor),
ExposedArea_shiftc = ExposedArea - c_ExposedArea * (Age - Age.cor),
K_shiftc = log10(TotalArea_shiftc) + 1/4*log10(AvgThickness_shiftc^2) - 5/4*log10(ExposedArea_shiftc),
I_shiftc = log10(TotalArea_shiftc) + log10(ExposedArea_shiftc) + log10(AvgThickness_shiftc^2),
S_shiftc = 3/2*log10(TotalArea_shiftc) + 3/4*log10(ExposedArea_shiftc) - 9/4*log10(AvgThickness_shiftc^2),
Knorm_shiftc = K_age_decay/sqrt(1 + (1/4)^2 + (5/4)^2),
Snorm_shiftc = S_age_decay/sqrt((3/2)^2 + (3/4)^2 + (9/4)^2),
Inorm_shiftc = I_age_decay/sqrt(1^2+1^2+1^2))
To compare the values, we will have the raw measurements in light gray, and the new values in dark blue
Comparing with the raw data:
Comparing with the raw data:
Combining the two processes is simple as extracting the slope from the linear mixed model used for the harmonization. Here, we also included the lobes.
# filter CTL subjects (again, if you are running this script with other diangostics)
data_rate <-
filter(data, Diagnostic == "CTL")
# assure Sample is a factor
data_rate$Sample <-
as.factor(data_rate$Sample)
# assure ROI is a factor and is in the desired order
data_rate$ROI <-
factor(data_rate$ROI,
levels = c("hemisphere", "F", "O", "P", "T"))
First, we run the lmer, then extract the residual for each sample and the *_shift* constant. After, we extract the Age.trend constant (called c in the previous deaging example). Finally, the code re-estimates the variable in all configurations: harmonized, deaged, and both.
Here, we used the linear version of the variables. It does not affect the final result.
## AvgThickness ---
m.1 <- lme4::lmer(AvgThickness ~ Age * ROI + (1 | Sample:ROI), data = data_rate)
re <- as_tibble(ranef(m.1)) %>%
filter(grpvar == "Sample:ROI") %>%
mutate(
T_shift = condval,
sd_T_shift = condsd,
Sample = str_split(grp,
pattern = ":", simplify = TRUE)[, 1],
ROI = str_split(grp,
pattern = ":", simplify = TRUE)[, 2]
) %>%
dplyr::select(-c(condval, grpvar, term, condsd, grp))
Age.trend <- as_tibble(lstrends(m.1, ~ ROI, var = "Age")) %>%
mutate(Age.trend_T = Age.trend) %>%
dplyr::select(c(ROI, Age.trend_T))
data <- full_join(data, Age.trend) %>%
full_join(re) %>%
mutate(
AvgThickness_shiftc = AvgThickness - T_shift,
AvgThickness_age_decay = AvgThickness - Age.trend_T * (Age - Age.cor),
AvgThickness_age_decay_shiftc = AvgThickness - T_shift - Age.trend_T * (Age - Age.cor)
)
## TotalArea ---
m.1 <- lme4::lmer(TotalArea ~ Age * ROI + (1 | Sample:ROI), data = data_rate)
re <- as_tibble(ranef(m.1)) %>%
filter(grpvar == "Sample:ROI") %>%
mutate(
AT_shift = condval,
sd_AT_shift = condsd,
Sample = str_split(grp,
pattern = ":", simplify = TRUE)[, 1],
ROI = str_split(grp,
pattern = ":", simplify = TRUE)[, 2]
) %>%
dplyr::select(-c(condval, grpvar, term, condsd, grp))
Age.trend <- as_tibble(lstrends(m.1, ~ ROI, var = "Age")) %>%
mutate(Age.trend_AT = Age.trend) %>%
dplyr::select(c(ROI, Age.trend_AT))
data <- full_join(data, re) %>%
full_join(Age.trend) %>%
mutate(
TotalArea_shiftc = TotalArea - AT_shift,
TotalArea_age_decay = TotalArea - Age.trend_AT * (Age - Age.cor),
TotalArea_age_decay_shiftc = TotalArea - AT_shift - Age.trend_AT * (Age - Age.cor))
## ExposedArea ---
m.1 <- lme4::lmer(ExposedArea ~ Age * ROI + (1 | Sample:ROI), data = data_rate)
re <- as_tibble(ranef(m.1)) %>%
filter(grpvar == "Sample:ROI") %>%
mutate(
AE_shift = condval,
sd_AE_shift = condsd,
Sample = str_split(grp,
pattern = ":", simplify = TRUE)[, 1],
ROI = str_split(grp,
pattern = ":", simplify = TRUE)[, 2]
) %>%
dplyr::select(-c(condval, grpvar, term, condsd, grp))
Age.trend <- as_tibble(lstrends(m.1, ~ ROI, var = "Age")) %>%
mutate(Age.trend_AE = Age.trend) %>%
dplyr::select(c(ROI, Age.trend_AE))
data <- full_join(data, re) %>%
full_join(Age.trend) %>%
mutate(
ExposedArea_shiftc = ExposedArea - AE_shift,
ExposedArea_age_decay = ExposedArea - Age.trend_AE * (Age - Age.cor),
ExposedArea_age_decay_shiftc = ExposedArea - AE_shift - Age.trend_AE * (Age - Age.cor)
)
Re-estimating K, S, and I:
data <- data %>%
mutate(
K_age_decay = log10(TotalArea_age_decay) + 1/4*log10(AvgThickness_age_decay^2) - 5/4*log10(ExposedArea_age_decay),
K_shiftc = log10(TotalArea_shiftc) + 1/4*log10(AvgThickness_shiftc^2) - 5/4*log10(ExposedArea_shiftc),
K_age_decay_shiftc = log10(TotalArea_age_decay_shiftc) + 1/4*log10(AvgThickness_age_decay_shiftc^2) - 5/4*log10(ExposedArea_age_decay_shiftc),
I_age_decay = log10(TotalArea_age_decay) + log10(ExposedArea_age_decay) + log10(AvgThickness_age_decay^2),
I_shiftc = log10(TotalArea_shiftc) + log10(ExposedArea_shiftc) + log10(AvgThickness_shiftc^2),
I_age_decay_shiftc = log10(TotalArea_age_decay_shiftc) + log10(ExposedArea_age_decay_shiftc) + log10(AvgThickness_age_decay_shiftc^2),
S_age_decay = 3/2*log10(TotalArea_age_decay) + 3/4*log10(ExposedArea_age_decay) - 9/4*log10(AvgThickness_age_decay^2),
S_shiftc = 3/2*log10(TotalArea_shiftc) + 3/4*log10(ExposedArea_shiftc) - 9/4*log10(AvgThickness_shiftc^2),
S_age_decay_shiftc = 3/2*log10(TotalArea_age_decay_shiftc) + 3/4*log10(ExposedArea_age_decay_shiftc) - 9/4*log10(AvgThickness_age_decay_shiftc^2),
Knorm_shiftc = K_shiftc / sqrt(1 + (1 / 4) ^ 2 + (5 / 2) ^ 2),
Snorm_shiftc = S_shiftc / sqrt((3 / 2) ^ 2 + (3 / 4) ^ 2 + (9 / 4) ^ 2),
Inorm_shiftc = I_shiftc / sqrt(1 ^ 2 + 1 ^ 2 + 1 ^ 2),
Knorm_age_decay = K_age_decay / sqrt(1 + (1 / 4) ^ 2 + (5 / 2) ^ 2),
Snorm_age_decay = S_age_decay / sqrt((3 / 2) ^ 2 + (3 / 4) ^ 2 + (9 / 4) ^ 2),
Inorm_age_decay = I_age_decay / sqrt(1 ^ 2 + 1 ^ 2 + 1 ^ 2),
Knorm_age_decay_shiftc = K_age_decay_shiftc / sqrt(1 + (1 / 4) ^ 2 + (5 / 2) ^ 2),
Snorm_age_decay_shiftc = S_age_decay_shiftc / sqrt((3 / 2) ^ 2 + (3 / 4) ^ 2 + (9 / 4) ^ 2),
Inorm_age_decay_shiftc = I_age_decay_shiftc / sqrt(1 ^ 2 + 1 ^ 2 + 1 ^ 2)
)
To compare the values, we will have the raw measurements in light gray, and the new values in dark blue.
ggplot(filter(data, ROI == "hemisphere")) +
geom_point(aes(x = Age, y = AvgThickness, alpha = 0.4), color = "grey") +
geom_point(aes(x = Age, y = AvgThickness_age_decay_shiftc, alpha = 0.4), color = "darkblue") +
geom_smooth(method = "lm", aes(x = Age, y = AvgThickness), color = "grey") +
geom_smooth(method = "lm", aes(x = Age, y = AvgThickness_age_decay_shiftc), color = "darkblue") +
guides(alpha= "none") +
theme_pubr()
ggplot(filter(data, ROI == "hemisphere")) +
geom_point(aes(x = Age, y = TotalArea, alpha = 0.4), color = "grey") +
geom_point(aes(x = Age, y = TotalArea_age_decay_shiftc, alpha = 0.4), color = "darkblue") +
geom_smooth(method = "lm", aes(x = Age, y = TotalArea), color = "grey") +
geom_smooth(method = "lm", aes(x = Age, y = TotalArea_age_decay_shiftc), color = "darkblue") +
guides(alpha= "none") +
theme_pubr()
ggplot(filter(data, ROI == "hemisphere")) +
geom_point(aes(x = Age, y = ExposedArea, alpha = 0.4), color = "grey") +
geom_point(aes(x = Age, y = ExposedArea_age_decay_shiftc, alpha = 0.4), color = "darkblue") +
geom_smooth(method = "lm", aes(x = Age, y = ExposedArea), color = "grey") +
geom_smooth(method = "lm", aes(x = Age, y = ExposedArea_age_decay_shiftc), color = "darkblue") +
guides(alpha= "none") +
theme_pubr()
t_a <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = AvgThickness, color = Sample, fill = Sample), alpha = 0.4) +
geom_point(aes(alpha =0.4)) +
geom_smooth(method = "lm") +
guides(alpha = "none") + # added some transparency to improve visualization
theme_pubr()+
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25))) # human lifespan
t_b <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = AvgThickness), alpha = 0.4) +
geom_smooth(method = "lm", fullrange = TRUE, aes( color = Sample, fill = Sample)) +
geom_smooth(method = "lm", fullrange = TRUE, color = "black", linetype = "dashed") +
guides(alpha = "none") + # added some transparency to improve visualization
theme_pubr() +
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25))) # human lifespan
t_c <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = AvgThickness_age_decay_shiftc, color = Sample, fill = Sample), alpha = 0.4) +
geom_point(aes(alpha =0.4)) +
geom_smooth(method = "lm") +
guides(alpha = "none") + # added some transparency to improve visualization
theme_pubr() +
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25))) # human lifespan
t_d <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = AvgThickness_age_decay_shiftc), alpha = 0.4) +
geom_smooth(method = "lm", fullrange = TRUE, aes( color = Sample, fill = Sample)) +
geom_smooth(method = "lm", fullrange = TRUE, color = "black", linetype = "dashed") +
guides(alpha = "none") + # added some transparency to improve visualization
theme_pubr() +
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25))) # human lifespan
ggarrange(t_a, t_b, t_c, t_d, common.legend = TRUE, ncol = 2)
## $`1`
##
## $`2`
##
## attr(,"class")
## [1] "list" "ggarrange"
ggplot(filter(data, ROI == "hemisphere")) +
geom_point(aes(x = Age, y = K, alpha = 0.4), color = "grey") +
geom_point(aes(x = Age, y = K_age_decay_shiftc, alpha = 0.4), color = "darkblue") +
geom_smooth(method = "lm", aes(x = Age, y = K), color = "grey") +
geom_smooth(method = "lm", aes(x = Age, y = K_age_decay_shiftc), color = "darkblue") +
guides(alpha= "none") +
theme_pubr()
ggplot(filter(data, ROI == "hemisphere")) +
geom_point(aes(x = Age, y = S, alpha = 0.4), color = "grey") +
geom_point(aes(x = Age, y = S_age_decay_shiftc, alpha = 0.4), color = "darkblue") +
geom_smooth(method = "lm", aes(x = Age, y = S), color = "grey") +
geom_smooth(method = "lm", aes(x = Age, y = S_age_decay_shiftc), color = "darkblue") +
guides(alpha= "none") +
theme_pubr()
ggplot(filter(data, ROI == "hemisphere")) +
geom_point(aes(x = Age, y = I, alpha = 0.4), color = "grey") +
geom_point(aes(x = Age, y = I_age_decay_shiftc, alpha = 0.4), color = "darkblue") +
geom_smooth(method = "lm", aes(x = Age, y = I), color = "grey") +
geom_smooth(method = "lm", aes(x = Age, y = I_age_decay_shiftc), color = "darkblue") +
guides(alpha= "none") +
theme_pubr()
k_a <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = K, color = Sample, fill = Sample), alpha = 0.4) +
geom_point(aes(alpha =0.4),
size = 4,
pch = 21,
colour = 'white') +
geom_smooth(method = "lm") +
guides(alpha = "none") + # added some transparency to improve visualization
theme_pubr() +
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25))) # human lifespan
k_b <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = K), alpha = 0.4) +
geom_smooth(method = "lm", fullrange = TRUE, aes( color = Sample, fill = Sample)) +
geom_smooth(method = "lm", fullrange = TRUE, color = "black", linetype = "dashed") +
guides(alpha = "none") + # added some transparency to improve visualization
theme_pubr()+
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25))) # human lifespan
k_c <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = K_age_decay_shiftc, color = Sample, fill = Sample), alpha = 0.4) +
geom_point(aes(alpha =0.4),
size = 4,
pch = 21,
colour = 'white') +
geom_smooth(method = "lm") +
guides(alpha = "none") + # added some transparency to improve visualization
theme_pubr() +
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25))) # human lifespan
k_d <- ggplot(filter(data, ROI == "hemisphere"), aes(x = Age, y = K_age_decay_shiftc), alpha = 0.4) +
geom_smooth(method = "lm", fullrange = TRUE, aes( color = Sample, fill = Sample)) +
geom_smooth(method = "lm", fullrange = TRUE, color = "black", linetype = "dashed") +
guides(alpha = "none") + # added some transparency to improve visualization
theme_pubr() +
scale_x_continuous(name ="Age",
limits=c(0,100),
breaks=c(seq(0,100,25))) # human lifespan
ap8 <- ggarrange(k_a, k_b, k_c, k_d, labels = c("A","B","C","D"), nrow =2, ncol = 2, font.label = list(size = 11), common.legend = TRUE, legend = "top")
ap8
ggsave("ap8.pdf", ap8, dpi=1200, width = 17.8, height = 17.8, units = "cm", device = "pdf")
This is by necessity a simplistic, one-parameter way of modeling aging; but this very simplicity allows a direct generalization to a broad set of morphological variables. In what follows we discuss this method’s known limitations.
Both processes are still dependable on a Control subset of the data. We use the control subset as the reference since the deaging, and the harmonization can remove unwanted effects if applied directly to nontypical or pathological populations. Although, it could be helpful to measure the difference in morphological variables due to different methodologies in the same Sample (e.g., comparing acquisition protocols).
We could use higher-level functions or split the data into periods of years to avoid considering a constant rate of aging for all human lifespan. Here we are limited by the available and processed data. International consortium and open datasets should fulfill this task. We also want to keep the processes as simple as possible.
We also expect future works to provide rates from smaller brain regions, leading to studies in precise ROI.